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In the present article we investigate optical near fields in semiconductor lasers. We perform finite element 
simulations for two different laser types, namely a super large optical waveguide (SLOW) laser, which is an 
edge emitter, and a vertical cavity surface emitting laser (VCSEL). We give the mathematical formulation of 
the different eigenvalue problems that arise for our examples and explain their numerical solution with the 
finite element method. Thereby, we also comment on the usage of transparent boundary conditions, which 
have to be applied to respect the exterior environment, e.g., the very large substrate and surrounding air. 

For the SLOW laser we compare the computed near fields to experimental data for different design pa- 
rameters of the device. For the VCSEL example a comparison to simplified ID mode calculations is carried 
out. 

Copyright line will be provided by the publisher 



1 Introduction The performance of novel nano-optical devices like SLOW lasers, or VCSELs depends 
very sensitively on the design of the system. Due to lack of practical experience and high fabrication costs 
of prototypes, numerical design and simulations are vital for the development of high performance devices. 
Numerical simulations can lead to optimized structures with, e.g., desired properties of the far or near field 
and also give insight into the physics of a device. 

Since in modern devices and applications the wavelength of light is of the same order as the dimension of 
the simulated structures, Maxwell's equations have to be solved rigorously, in order to get accurate results 
for the electromagnetic field. Examples of such structures are semiconductor lasers, meta materials, pho- 
tonic crystal devices, photolithographic masks and nano-resonators. A lot of different simulation techniques 
have been applied to and developed for nano-optical simulations, e.g., the finite element method (FEM), 
finite difference time domain simulations (FDTD), wavelet methods, finite integration technique (FIT), or 
rigorously coupled wave analysis (RCWA) |4|. 

Here we present the finite element method for the solution of time harmonic Maxwell's equations in 
semiconductor laser devices. The finite element method offers the possibility for high order discretization 
schemes, which leads to very accurate results at low computational costs. In this work the FEM solver 
JCMsuite developed at the Zuse-Institute Berlin and JCMwave GmbH was used for simulations l9l [TT1 l6l . 
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(b) 




Fig. 1 (a) Geometry of a propagating mode problem with invariance in z-direction; (b) corresponding 2D cross section. 
The geometry is divided into an interior domain Qi n with boundary T. The structured exterior domain extends to infinity. 



2 Mathematical formulation of eigenmode problems In our numerical analysis of semiconductor lasers 
we focus on an edge emitting laser and a VCSEL. Our goal is to compute optical near field distributions of 
the lasing modes within these devices. In practice these simulations can be used for the design of real world 
devices, e.g., to study the dependence of the light pattern emitted from the laser on geometrical parameters. 
As an illustration such an analysis will be carried out for an edge emitting laser and we will compare our 
results to experimental measurements in Section 15.21 The computed lasing eigenmodes are also a crucial 
ingredient for full electro-optical device simulation |3). However, this will be part of future work. 

Mathematically computation of modes for the edge emitting laser and the VCSEL example lead to differ- 
ent eigenvalue problems. For the SLOW laser a propagating mode problem has to be solved and a resonance 
problem for the VCSEL. We are not only interested in the near field distribution, but also want to determine 
leakage losses of the lasing modes. Therefore, we can not restrict the domain of interest, e.g., to a cavity of 
finite size. The eigenvalue problems are per se stated on unbounded domains, which take into account the in- 
finite surrounding exterior, i.e., the substrate and surrounding air. Since in a numerical simulation the size of 
the computational domain has to be restricted, we have to apply appropriate transparent/radiating boundary 
conditions. This is done by formulating the eigenvalue problem as a coupled interior-exterior problem. 

We are interested in stationary solutions for the electric field E. Therefore, we make the following 
harmonic ansatz for the time dependence: 



E(x,y,z,t) = e- iut E(x,y,z), 



(1) 



where tu denotes the frequency. Using this ansatz in Maxwell's equations, the following second order "curl- 
curl equation" for the electric field can be derived: 



x E-u/eE = 0, 



(2) 



where we assumed no exterior sources. The quantities /i and e denote the permeability and permittivity 
tensor respectively. 

This formulation is used in the following for derivation of the eigenvalue problems. 



2.1 Propagating mode problem Propagating mode problems arise, when waveguide structures as de- 
picted in Fig. Q2a) are considered. The geometry has an invariance in one spatial dimension (here we choose 
the z-direction). Then the following ansatz for the electric field can be made: 

E(x,y,z)=e ik * z E(x,y), (3) 
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where k z is the propagation constant. Hence, we assume a harmonic dependence of the electric field in the 
invariant direction. Note that with ansatz (|3), the differential operator with respect to z can be replaced as 
follows: 



d z = ik z . (4) 

With ansatz ([3), furthermore, we only have to consider the cross section of the waveguide in the xy-plane. 
Figure Q2b) depicts this infinite domain of interest. The domain is divided into a finite interior S!j n C 2 and 
an infinite exterior 2 \ £!;„■ The boundary of Sl; n is denoted by V. The propagating mode problem can now 
be given as follows: 

Problem 1 (Propagating mode problem) For given frequency u find pairs ([Ei n , E cx ] , k z ) such that: 
(i) The electric field E; n fulfills Maxwell's equations in the interior domain: 



d x s 




9y 




ik z 


1 \ Cti z 



X 



E in (a;, y) - cj 2 eE in (x, y) = in O in . (5) 



(ii) The electric field E cx fulfills Maxwell's equations in the exterior domain: 

x E ex (a;, y) - w 2 eE cx (x, y) = in 2 \ Cl in . (6) 



d x N 




' d x 


8y 




8y 


ik z , 




. ik z 



(iii) Boundary condition at T: the tangential component of the electric field is continuous: 

n x (E in - E ox )| r = 0, 
where n denotes the outward pointing normal vector on T. 

(iv) Radiating boundary condition: E cx is strictly outward radiating. 

For homogeneous exteriors this can be stated, e.g., by the Silver-Miiller radiation condition: 

lim r ( V x E sc (r) x r - ^V e ° xt ^ cxt y x E sc (r) ) = 

i — >cx> \ C J 

uniformly continuous in each direction ro, 

where r is the coordinate vector in 2 , r its norm, r = -, and e ext and /i ext the relative permittivity and 
permeability in the exterior. 

A more general characterization of outward radiating fields for structured exterior domains is the pole 
condition lfT0ll5l. 



Problem Q] yields a quadratic eigenvalue problem for the electric field and the propagation constant k z . The 
division of the electric field into an interior and exterior field has to be done, in order to guarantee that 
the computed modes are outward radiating. Physically this means that electromagnetic field energy of a 
mode within the waveguide only propagates outwards. Since the propagating mode problem is stated on an 
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unbounded domain, the eigenvalue k z becomes complex. Hence, the propagating mode looses energy while 
propagating in z-direction and is damped according to: 

E(x,y,z) = e -^ fe *>* (e m{k * }z -E(x,y)) , 



where di{k z } and denotes the real and imaginary part of the propagation constant k z . In numerical 

simulations we will compute the effective refractive index: 

"oft = 7^, (7) 
instead of the propagation constant k z . Thereby k = oj/c is the vacuum wave number. 



2.2 Resonance problem The setup of a resonance problem differs from the propagating mode problem. 
Here, the electric field and its resonance frequency have to be determined for a given geometry, which can 
be 2D or 3D. Again the domain of interest is unbounded, and it is divided into an interior 0; n C 3 and an 
exterior ^i n \ 3 , analogously to Fig. Q] The coupled interior-exterior resonance problem is given as follows: 

Problem 2 (Resonance mode problem) Find pairs ([Ei n , E ex ] , uj) such that: 

(i) The electric field Ej n fulfills Maxwell's equations in the interior domain: 

V x x E ia (x, y, z) - o; 2 eE in (a;, y, z) = in Q in . (8) 

(ii) The electric field E cx fulfills Maxwell's equations in the exterior domain: 

V x fi- 1 V x E ex (x, y, z) - w 2 eE cx (x, y, z) = in 3 \ O in . (9) 



(iii) Boundary condition at T: the tangential component of the electric field is continuous: 

nx (E in -E cx )| r = 0. 



(iv) Radiating boundary condition: E cx is strictly outward radiating. 



Problem |2] is an eigenvalue problem for the frequency ui 2 and the electric field. Again the division of the 
electric field is necessary to obtain outward radiating resonances. Since the resonance problem is stated on 
an infinite domain, there are no real-valued eigenvalues. The eigenfrequencies become complex, and their 
imaginary parts correspond to the inverse lifetimes and radiation losses of the corresponding eigenmodes: 

E(x, y, z, t) = e~ Tt (e^'E^, y, zj) , with r = -— 



3 Transparent boundary conditions The eigenvalue problems introduced in the previous section are 
stated on unbounded domains. This makes straightforward discretization impossible, since we can only 
consider finite domains in a numerical simulation. Here we use the perfectly matched layer (PML) method 
(2) to reduce the domain of interest to finite size. Hence, the PML can be seen as a transparent boundary 
condition. 
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Fig. 2 First order vectorial ansatz functions defined on triangle of finite element discretization. 



.WW, 

,w\\\, 
■ w\\\\. 

,\\\\\\\, 
, \ \ \ \ \ \ \ } 

ililllii 



0.5 



The key idea of the PML method is a complex continuation of the electric field. Let us consider a ID 
setup, where the positive real line is the infinite exterior. An outward radiating field, which travels from the 
origin x = to the right, is given by a plane wave: 



E(x) 



with k > 0. 



Now we introduce a new coordinate system as follows: 

x = (1 + ia)x , with < a £ . 
The plane wave in the new coordinate system is given by: 

E(x) = e- CTfc V fc5 , 

hence, it is decaying exponentially with increasing distance from the origin. Therefore, we can set E(x) = 
for a sufficiently large x. This means we can truncate the exterior domain at a finite value. 

In 2D and 3D the implementation of this idea is more technical lfl2l . First a coordinate system in the 
exterior has to be introduced, which includes a generalized distance variable £. Then a complex continuation 
along this distance variable is carried out, and the exterior domain is truncated after sufficient decay of the 
complex continuation of the exterior field. At this artificial boundary the electric field or its derivative are 
set to zero, corresponding to zero Dirichlet or Neumann conditions. This guarantees that the exterior field is 
strictly outward radiating. 



4 Finite element discretization For finite element discretization the curl-curl equation for the electric 
field © has to be stated in so-called weak form: 
Find E e H (curl , 0), such that: 

o(E, v) = f(v) , VdgH (curl , Q) . 

This is a variational formulation of Maxwell's equations, where the electric field E is determined from a 
proper infinite dimensional function space H (curl , f2) Q, and the bilinear form a(-, •) basically represents 
the curl-curl equation. This offers a very elegant discretization scheme. The finite element version of 
Maxwell's equations simply reads: 
Find Eg 14, such that: 

o(E, ») = /(«), VveV h , 
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Fig. 3 Cross section of cylinder symmetrical VCSEL. 

where the variational formulation is stated on a finite dimensional function space Vh C H (curl , 0), with 
dim Vh = N < oo. The finite element space Vh is constructed as follows. First the domain of interest is 
subdivided into a number of small patches, i.e., it is triangulated. Then on each patch a basis of polynomial 
ansatz functions of maximum polynomial degree p is chosen, see Fig. [2] The union of all ansatz spaces on 
all patches gives the finite element space. 

The quality of a finite element solution can be increased by refining the patches, i.e., the mesh. This is 
referred to as /i-refinement, where h is the measure of the dimension of the patches. Another alternative is 
increasing the polynomial degree of the ansatz functions on each patch, which is referred to as p-refinement. 
For smooth solutions in general p-refinement (high order) leads to much faster convergence. We will analyze 
this in numerical examples. 

5 Finite element simulation of semiconductor lasers In this section we use the finite element method 
for computation of lasing modes of two types of semiconductor lasers. The first example is a 3D cylinder 
symmetrical VCSEL. We will compare our results to ID simulations, which are often used in the design 
process due to their computational simplicity. 

The second example is a SLOW laser, which is produced at the Ferdinand Braun Institute (FBH). For 
this device experimental data of the near field distributions of the lasing modes is available, which can 
be compared to numerical results. The influence of design parameters is studied, and the finite element 
simulations will provide an explanation of the experimentally observed behaviour. Further we will analyze 
the convergence of the finite element simulation for this example. 

For the waveguide problem Q] the frequency lu of the mode is fixed. Hence the dispersion of materials 
is not important. For the resonance mode problem the frequency lu is the unknown eigenvalue. Since the 
material parameters e depend on to this leads to a non-linear eigenvalue problem in general, which can 
be solved iteratively. Usually the DBR mirrors of a VCSEL are designed for a specific frequency. This 
frequency is used as an initial guess to fix the material parameters. After the eigenvalue lu is computed, 
the material parameters are determined for the new frequency value. This procedure is continued until the 
eigenvalue lu is converged. 

5.1 Vertical cavity surface emitting laser The layout of the VCSEL example is shown in Fig. [3] The 
InGaAs active layer is embedded into a top and bottom GaAs/AlGaAs DBR mirror. Two AlOx appertures are 
placed above the active zone. The corresponding material parameters were computed from [l] for AlGaAs 
and taken from Landolt-Bornstein for AlOx and InGaAs. 
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Fig. 4 Intensity of electric field obtained from FEM simulation in linear (top) and logarithmic (bottom) scale. In 
logarithmic scale the out-coupled field is visible. 

Since the design is cylinder symmetrical, simulations can be carried out on a 2D cross section. Figure 
|4] shows the fundamental mode of the VCSEL. On logarithmic scale we observe that the resonating mode 
radiates to the exterior. The corresponding resonance frequency is given by: 

w 3D = 1.926 • 10 1S + 1.86 • 10 n i, 

which corresponds to a wavelength of: 

A 3D = 978.12 nm. 

The Bragg mirror was designed for a wavelength of 980 nm. 

For the ID computations a scattering problem for the dielectric profile along the VCSEL's symmetry 
axis was solved with varying wavelength of an incidence plane wave. The resonance wavelength was then 
determined from a reflectivity dip of the ID profile, see Fig. |3a). We found: 

Aid = 979.12 nm, 

which is lnm apart from the 3D value. For comparison of the 3D model to ID calculations, we extract 
the finite element solution along the symmetry axis of the computational domain. Figure [2b) shows the 
corresponding permittivity profile along this axis together with the electric field intensity from ID and 3D 
simulations. The near fields agree to a large extend. Of course the ID simulations do not give the lateral 
mode shape and can not be used for lateral design of the VCSEL. 

5.2 Super large optical waveguide laser The SLOW laser fabricated at the FBH was designed for a wave- 
length of A = 1060 nm, which was used in numerical computations. It consists of a 9 /im broad vertical 
waveguide core in the middle of which the active region is placed and a wide ridge waveguide defined by 
lateral trenches filled with electroplated gold. The relative permittivities of the waveguide used in our simu- 
lations was e = 11.7. The permittivities for the top and bottom claddings was fixed at e = 11.3. More details 
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Fig. 5 (a) Reflectivity of ID VCSEL model in dependence on incidence wavelength. The permittivity profile is given 
in (b): permittivity profile and electric field intensity of the fundamental mode along symmetry the axis of the VCSEL. 
A comparison of ID semi-analytic and 3D finite element simulation is shown. 



Fig. 6 Finite element triangulation of cross section of super large optical waveguide laser. The meaning of the parame- 
ters w and d is also indicated. 



on the design and experimental results can be found in J8). The widths of both the vertical and the lateral 
waveguides are by a factor of nearly 10 larger than is commonly used. The design parameters are, e.g., the 
width w of the ridge and the residual layer thickness d, which denotes the distance between each trench and 
the active layer. Etching the trenches deeper into the material, reduces the residual layer thickness d. The 
finite element triangulation of our second example is shown in Fig. [6] 

Let us look at the accuracy and computational costs for this example. Figure shows the convergence of 
the effective refractive index for different finite element orders p in dependence on the spatial discretization. 
The convergence for all orders p was computed against a reference value of: 

n eS = 3.422595157+ 1.82097- 10~ 4 i. 

This value was obtained from the most accurate finite element solution with order p = 4 and a spatial 
discretization leading to N = 3956993 finite element degrees of freedom. We observe that for high order p, 
less finite element degrees of freedom are necessary, to obtain an accurate solution. However, for all orders 
the relative error of the eigenvalue is very small, below 10~ 4 and down to 2 ■ 10 -8 for finest discretization. 
For order p = 1 first the relative error rises, and the convergent regime starts at N ~ 30000. Figure [7] also 
gives single CPU times, where 10 modes closest to an initial guess were computed. Within a few minutes, 
very accurate finite element solutions can be computed for this real-world example. 
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Fig. 7 Convergence of fundamental eigenvalue for SLOW laser example. The relative error against a very accurate 
solution (computed with higher spatial discretization and p = 4) is shown in dependence on number of finite element 
degrees of freedom (uniform refinement of triangulation) for different finite element orders p. Single CPU times are 
given for orders 2 and 4 (the computations were performed on Dual-Core AMD Opteron Processors with 2.8 GHz) 
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Fig. 8 (a) Experimental and (b) numerical lateral near field "intensity profiles" of fundamental lasing mode for different 
values of w, given in fim. 



Figure[8ja) shows corresponding experimental lateral near field profiles for different values of the design 
parameter w of the SLOW laser. The finite element results for the same design parameters are given in 
Fig. [HJb), where the data are obtained from a lateral cut through the active zone. In order to compare 
experimental and numerical data, we extract the widths (FWHM, full width at half maximum) of the near 
field intensity profiles in dependence on the parameter w. The result is depicted in Fig. |9|a). We observe a 
distinct minimum of the near field width for a value of w = 5 fjm. Although the experimental and numerical 
data show an offset, which is probably caused by the resolution limit of the measurement setup, the position 
of the minimum agrees. We observe that the lateral dimension of the lasing mode increases with increasing 
width of the ridge. However, reducing the width below w = 5 /im leads to a broader near field distribution. 
In order to explain this behaviour, the near field distributions obtained from finite element computation for 
different values of w are given in Fig. [10] We observe that for large w the lasing mode basically fills the 
ridge. For very small width the lasing mode, however, is no longer confined by the ridge and extends below 
the trenches towards the outer region. 
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Fig. 9 (a) Comparison of the widths (FWHM) of the experimental and numerical lateral near field intensity profiles 
corresponding to Fig. [8j (b) widths of the numerical lateral near field intensity profiles in dependence on residual layer 
thickness d for different values of w. 



In Figure |9jb) we analyze the dependence of the lateral width of the near field on the residual layer 
thickness d. We observe that for large w the residual layer thickness only has a weak influence on the width 
of the near field. However, for smaller values of w, a decreasing residual layer thickness leads to stronger 
lateral confinement of the lasing mode. Again finite element simulations help to understand this behaviour. 
Since for large w the lasing mode is mainly located between the trenches, c.f. Fig. [10] the residual layer 
thickness does not effect the mode. For a small value of w the mode extends below the trenches, and a 
decreasing residual layer thickness basically leaves no space for the mode below the trenches. The mode is 
then again confined between the trenches, see Fig. QT| 

6 Conclusions We showed that the finite element method is very well suited for the computation of 
optical modes in semiconductor lasers. A convergence analysis for an edge emitting laser demonstrated that 
accurate numerical results can be obtained at low computational costs. 

The numerical results were compared to real-world lasers. The dependence of experimentally obtained 
near field distributions on design parameters of the lasers could be reproduced and understood with the help 
of simulations. 

Furthermore, the FEM results for a 3D VCSEL were compared to a simplified ID model. The near field 
distributions of both models, thereby, agreed very well, whereas the computed eigenvalue differed slightly. 

7 Acknowledgment This work was carried out within SFB787 "Halbleiter Nanophotonik", funded by 
the DFG. 
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1.0, 2.0/im. Gold trenches and active layer are shown in green and red respectively. 
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